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Abstract 

The spectral test of random number generators (R.R. Coveyou and R.D. McPher- 
son, 1967) is generalized. The sequence of random numbers is analyzed explicitly, 
not just via their n-tupel distributions. The generalized analysis of many generators 
becomes possible due to a theorem on the harmonic analysis of multiplicative groups 
of residue class rings. We find that the mixed multiplicative generator with power 
of two modulus does not pass the extended test with an ideal result. Best qualities 
has a new generator with the recursion formula X^^i = aXk + cint(/c/2) mod 2°^. 
We discuss the choice of the parameters a, c for very large moduli 2'^ and present an 
implementation of the suggested generator with d = 256, a = 2^^^+2^^+2'^^+62181, 
c = (2160 + 1) . 11463. 
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1 Introduction 

The spectral test was proposed by R.R. Coveyou and R.D. McPherson in 1967 [|T|. The 
advantage of this test is to present an algebraic criterion for the quality of the generator. 
For the mixed multiplicative generator X^+i = aXk + cmodM 

min{|s| = \J Si + . . . + s"^ with Sa = Si + as2 + ■ ■ ■ + a!^~^Sn = OmodM} (1) 

should be as large as possible |l|, 0. The criterion is that simple since the n-tupels of 
random numbers form an n-dimensional lattice (cf. e.g. Fig. 4). A good generator has 
uniformly distributed ra-tupels which refers to an almost cubic lattice 0, ^, ^. 

The lattice is a consequence of the (affine) linear dependence of Xk+i on X^. From 
the figures on the left (type I) we see that the relation between k and X^ is much more 
complicated. This is however one of the most fundamental aspects of randomness. In order 
to judge whether a sequence X^ takes random values one would first plot the sequence 
itself and then maybe Xk+i over X^. 

Of course, the correlation between k and X^ is not independent from the distribution 
of pairs (X^, X^+i). E.g., a poor 'random' sequence X^ = ak lying on a line with gradient 
a leads to pairs (Xfc,Xfc+i = X^ + a) lying on a fine with gradient 1 shifted by a off 
the origin. This makes it reasonable to judge randomness by only looking at the n-tupel 
distributions. However, random number generators which have identical valuation by the 
spectral test may still look quite different. The generators X^+i = 41Xfc + 3 mod 1024 
(Fig. 5) and X^+i = 41Xfc + 1 mod 1024 (Fig. 6), e.g., differ only by the additive constant 
which does not enter Eq. (|I]). The lattices of pair distributions (type II in the figures) 
are similar whereas the plots of X^ over k show different behavior. The spectral test not 
even makes a difference between a prime number and a power of two modulus (cf. Fig. 1 
vs. Fig. 4). 

Therefore it is desirable to include the analysis of the correlation between k and Xk 
into the valuation of the test. In fact it is possible to analyze the accumulation of random 
numbers along certain lines (which is often seen in the figures) by Fourier transformation. 
More generally we extend the spectral test by analyzing the correlation between k and 
the n-tupel (Xfc,Xfc+i, . . .,Xfc+„-i). The generators mentioned above (Figs. 5, 6) acquire 
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different valuations. Fig. 6 is preferred since the random numbers spread more uniformly 
in Fig. 61 than in Fig. 51 (cf. Ex. 4.3 l.(a)). 

We will develop a powerful theorem (Thm. 3.4, Cor. 3.5) on the harmonic analysis 
of multiplicative groups of residue class rings that allows us to algebraically perform the 
test for the standard generators. It becomes even possible to evaluate the extended test 
on more generators than the original spectral test has yet been applied. This has the 
advantage to have more freedom for the search after an optimum generator. 

In fact we will find that the commonly used mixed multiplicative generator always 
shows correlations along certain lines if the modulus is not a prime number. We will 
improve the generator until these correlations will essentially disappear (from Fig. 4 via 
Fig. 9 to Fig. 11). We finally recommend the recursion formula 

Xk+i = aXk + c int {k/2) mod 2'^ (2) 

with the parameters 

d = 2^ do , a = 2^"'"'^'' + 2^"'^'^'' + ... + 2^"^° + (3 580 621 541 mod2'^«) , 

c = (2^^"^'/^^^ + 1^ (3 370 134 727mod2'^«) . (3) 

In particular, the case do = 16, k = 4 is discussed in Ex. 5.1 1. 

This generator is supposed to be a good choice with respect to the following three 
criteria. 

Firstly, the sequence of numbers provided by the generator should behave as close to 
a true random sequence as possible. 

Secondly, the calculation of random numbers should be as fast as possible. The gen- 
erator given in Eqs. @, is explicitly constructed to have best performance. It is 
important to note that this is not independent from the first criterion. It is possible to 
produce better random numbers the more effort one spends in calculating the numbers. 
Figs. 4 and 8 show how a simple doubling of the digits of the modulus improves the ran- 
domness of the generator. In general we can produce arbitrarily good random numbers 
with e.g. do = 32 and large k in Eq. (|). 

Thirdly, the properties of the random numbers should be known as detailed as possible. 
It is not sufficient to use a messy, opaque formula. It has often been seen that this leads 
to numbers which are far from being random 0. As long as one is not familiar with 
the qualities of the generator one can never rely on the results gained with it. The full 
evaluation of the generalized spectral test is supposed to provide a profound knowledge 
of the generator. 

We start with the development of the generalized spectral test in the next section. In 
Sec. D we derive some mathematical results on the Fourier transformation of residue class 
rings. In Sec. ^ we apply the result of the previous section, Thm. 3.4, to a series of 
commonly used and some new generators. Finally we discuss the choice of parameters in 
Sec. §. 
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2 The generalized spectral test 
2.1 Review of the spectral test 

We start with a short review of the spectral test [|1], 0] in which we try to stress its 
geometrical meaning. The idea is to plot all n-tupels of successive random numbers in an 
ra-dimensional diagram. This is done, e.g. for n = 2 in the figures of type II. 

Mathematically a figure is presented as a function g which is 1 at every dot and 
elsewhere. If A'x is the period of the generator X, that is the smallest number with 
Xk+Nx = V/c, then 

Nx 

k=l fcgZjv^ 

where 6a,b = 1 if a = 6 and 6a,b = if a 7^ 6 (for later convenience we also write the 
Kronecker 6 as Sa=b)- Moreover we have introduced the notation 

Xfc = (Xfc,Xfc+i, . . .,Xfc+„_i) , X = (xi,X2, . . , Znx = '^/NxZ. (5) 

We want to check whether the dots accumulate along certain hyper-planes (see e.g. 
the lines in Fig. III). To this end we select a hyper-plane and project all the dots onto a 
line perpendicular to it. If points accumulate along the plane many dots will lie on top 
of each other, otherwise the dots are spread uniformly over the line. 

The hyper-plane H is determined by its Hesse normal form 



= {x : SiXi + S2X2 + . . . + SnXn = s ■ X = 0} , |s| = s ■ s = ^ sj + . . . + ^ . (6) 

The line is stretched by the factor |s|. The position of a point on the line perpendicular 
to the plane is given by the number s ■ X^. 

Next we wind the line up to a circle so that the modulus M as point on the line lies 
on top of the 0. For a suitable choice of s, namely s = (—3, 7) all points in Fig. Ill lie 
now on the point represented by the number 256. 

The points are realized as complex phases on the unit circle. We obtain the assignment 

Xfc exp [ — s ■ Xfcj . (7) 

Finally we draw arrows from the center of the circle to all the dots and add them. The 
length of the resulting vector describes how the dots are balanced on the circle. If the 
dots spread uniformly the arrows cancel each other and the resulting vector is small. If, 
on the other hand, all dots lie on top of each other the length of the arrows sums up to 
Nx. 

If we restrict ourselves to integer si, . . . , s„ (accumulation of random numbers always 
occur along hyper-planes given by integer Si) the resulting vector is given by the Fourier 
transform of g, 
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where we have introduced the normahzation factor Nx The information about accu- 
mulations along the hyper-plane is contained in the phase of g is irrelevant. 

We remember that for the mixed multiplicative generator the n-tupels form a lattice 
(which is displaced off the origin). So |^P(s) will assume the maximum value Nx if s 
lies in the dual lattice, s ■ = C + iM, C,i E Z, otherwise |5'P(s) is zero. Since 
Xk = c{a^ — l)/(a — 1) modM this means, if gcd(c, M) = 1 and X has full period, that 
VA; : {a'' - l)sa = OmodMgcd(a - 1, M) from which Eq. (|ip follows (cf. Ex. 4.3 1). 



2.2 Generalization of the spectral test 

We generalize the spectral test by caring for the sequence in which the n-tupels are 
generated. The index k is added to the n-tupel X^ as zeroth component and we define g 

as 

^(xo,x)= ^^o,fc^x,Xfe = 5x,x., • (9) 

The geometrical interpretation remains untouched but now we consider also the figures of 
type I. With the Fourier transform of g we introduce some notation that will be needed 
later, 

g^,^ [iV2, M, X] (.0, s) = E exp (— .0^ + ■ X, j . (10) 

Normally A^i = A^2 = Nx is the period of X. We write g[N2,M,X] if Ni = N2 and we 
also suppress A^2 if Ni = N2 = Nx- However, sometimes it is convenient to use a multiple 
N2 of the period Nx instead of the period itself. Moreover for some values of (sq, s), e.g. 
(0, 0), the right hand side is periodic in k with a smaller period Ni\N2. In this case we 
may sum over only and use the subscript I^Ni- Note that Z]fceZ]v /(^) always implies 
that / is a function on Zat^ which means that / is periodic, f{k) = f{k + A''i) Vfc. We 
also suppress M and X on the left hand if the context is clear. 

The sum over k is hard to evaluate since in the exponential k is combined with X^- 
However in fact we are interested in \g\'^ and find 

\9\Hso,s) = E exp (l^so (/c' - A;) + ■ (Xfc, - Xfc) 

fc.fc'ezjv^ ^Nx M 

The sum over k has no linear fc-dependence, only differences of random numbers occur. 
Like in the standard spectral test in many cases the sum over k can be evaluated. The 
result is often simple enough to be able to evaluate the sum over Ak also. 

If, e.g., Xk = a'' modP where P is a prime number and a is a primitive element of Zp, 
the multiplicative group of Zp, we find that s- (X^+aa; — X^.) = Saa'^(a^'^ — 1) = Sak{a^^ — 
1) modP for some 7^ G Zp. If k runs through the P — 1 values of Zp then k sweeps 
out the whole Zp\{0}. The sum over k can be evaluated yielding P6/\k=o — 1 for Sa 7^ 0. 
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Finally the sum over Ak gives together with the normalization = P/{P — 1) — 6so=o 
(cf. Eq. (|5|) for c/= 1). 

Note that the standard spectral test corresponds to sq = 0. We give some simple 
results on l^fp in the following lemma. 
Lemma 2.1. 

\g\^[ciXk+c2 + cs] (so, s) = I^HXfc] (so, cis) 

= I^HXfc] (so,s) if 3co : ciX^ = Xk+co , (12) 
\9\l,.^j,JciNx]{so,s) = C2(5,o=omodcil^P (so/ci,s) , (13) 

l^lli = 1 > Igl'^ iso,0) = Nx6so=OmodNx , (14) 

E |^r(so,s) = iVx , E l^l'(^o,s) = M"ifXfc = Xfc, ^fc = fc'modiV^. (15) 

Proof. The proofs are straight forward. We show Eq. (^) to get used to the notation: 

2-711 . 2'Ki 



^Z.^iV;, [CliVx] (So, S) = ^l/c2Nx E I ^r]\^^o/^ + ■ Xfe 



/ 9Ar / 27ri , 2ni 

C2I cfNx 2^ exp [ -^r^^ok + — s • Xfe 

ca/cfA^x E E exp (^-^so (fci + AT^^fca) + ■ Xfc, 



C2/ iVx5so=o mod ci E — ^1 + irr^ ■ -^^=1 



2Tli Sq 2TTi 



Nxci M 



□ 

One may also be interested in correlations between non-successive random numbers 
like Xfc and In general it is possible to study n-tupels ^k+r = {^k+n, ■ ■ -y^k+Tn)- 

This amounts to replacing X^ by ^k+r and Sa by Sa,T = Sia'^^ + . . . + s^a"^" in our results. 



2.3 Valuation with the generalized spectral test 

Now we have to clarify how the calculation of leads to a valuation of the generator. 

We can not expect that \g\'^ vanishes identically outside the origin since in this case g 
would be constant. Eq. (|1^) shows that the mean value of is 1. 

What would we expect for a sum of truly random phases? Real and imaginary part 
of a random arrow with length 1 have equal variance 1/2. For large Nx the sum of 
arrows is therefore normally distributed with density l/nNx ■exp(— (x^ + y'^)/Nx)dxdy = 
exp{—r'^/Nx)dr'^/Nx- Thus z = has the density exp(— z) for a true random sequence 
and the expected value for \g\'^ is 1. 

This means that values of < 1 can be accepted. It is clear that for a given 
(so,s) 7^ (0,0) the correlations are worse the higher \g\'^{so,s) > 1 is. But what does the 
location of an (sq, s) with |5fp(so, s) > 1 mean for the generator? 
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We remember that (sq, s) may be seen as normal vector on the hyper-plane along which 
the accumulations occur. If e.g. n = 1 and (sq, si) = (1, 1) the corresponding 1-plane has 
the equation xo + xi = (cf. e.g. Figs. 41, 51). If the fc-axis and the X^-axis are normalized 
to length 1 this line has length a/2. With the normal vector (3,1) (cf. Fig. 61) one obtains 
the equation 3xo + xi = mod 1 which intersects the unit cube three times and therefore 
has the length a/S^ + 1 = v^TO. Accumulations along this longer line are less important 
than along the short line. In the extreme case where the line fills the whole unit cube 
by intersecting it very often, accumulations can hardly be recognized. Note that in this 
sense the normal vectors (sq, si) and (2so, 2si) do not determine the same hne. The latter 
one contains e.g. the points (1/2,0), (0, 1/2), (1, 1/2), (1/2, 1). It has twice the length of 
the former one and too large a has half the effect. 

We generalize these considerations to n > 1 by taking the area of the n-dimensional 
hyper-plane with normal vector {sq, s) as measure for the importance of the accumulations 
detected. The area is given by |(so, s)| = (sq + s^)^/^, the Euclidean length of the normal 
vector. 

We can relate both mechanisms by defining the quality parameter 

I ('^0 s) I 

Qn (go, s) = ' , Q„ = max(^ s)ga xKA{o,o}Qn (so,s) . (16) 
\g [So, sj p ' X M J 

Good generators have Qi ~ 1. It is hard to achieve Qn ~ 1 for n > 1 (see however Sec. 
[4.5|) . More realistic is Qn ~ M^/"~^ (cf. Sec. |^) which means that the distribution of 
n-tupels deteriorates for higher n. In general small n are more important than large n. 
Apart from the value of Qn also the number of sites (so,s) at which Qn{so,s) = Qn is 
relevant (cf. Ex. 4.3 l.(a)). 

Let us try to find an interpretation for Q„. Assume the generator produces only mul- 
tiples of t\M. Then |^p(0, Si = M/t, 0, . . ., 0) = Nx, thus Q„(0, M/t, 0, . . ., 0) = M/tNx, 
and NxQn determines the number of non-trivial digits. In general NxQn may be larger 
than M and therefore we say that M„ = max(AOf Qn, M) determines the number of digits 
we can rely on. Analogously Nn = max(A'xQn, Nx) gives the quantity of random num- 
bers for which the n-tupel distributions are reasonably random. Specifically Mn(so,s) = 
max(A^x<5n('5o, s), M) determines the digits and iV'„(so,s) = max(A'xQn(so, s), A^x) the 
quantity of random numbers not affected by accumulations perpendicular to (sq, s) (cf. 
0, ,P- 90]). 

Note however that these are only crude statements. If, e.g., the 'period' of the gen- 
erator is enlarged by simply repeating it then NxQn{sQ,s) remains unaffected only if 
So = 0. Moreover, a high |^p(so,s) may be harmful even if |(so,s)| is large (cf. Fig. 31 
with Qi > 1). 

Note that Qn is a relative quality parameter. Although Qn usually does not increase 
with larger modulus (for > 1 is actually decreases) the quality of the generator gets 
better since Nx grows (cf. Fig. 4 vs. Fig. 8). 
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3 Multiplicative groups of residue class rings 

This section provides the mathematical tools that are needed to analyze several random 
number generators. We derive Eq. which determines for = a'^modM for 
various moduli M and multipliers a. There are some connections to the theory of Gaufi 
sums (see e.g. [P, §3]) where however the calculation of the absolute is almost trivial. In 
our case the exponent is not quadratic but has a linear and an exponential component 
(cf. Cor. 3.5, but also Ex. 4.6 4). The modulus is not restricted to a prime number. We 
are concerned with residue class rings rather than with fields. 
Throughout the paper we use the following notation, 

Zm = Z/^Af integers modM. 

Z]^^ = multiplicative group of , a E Z]^j <^==^ gcd (a, M) = 1 . 
{a)M = {a^modM, G Z} , the cychc group of a modM . 

Sa = E"=lSja-'"S Sa,M = gcd (Sa, M) . (17) 

N = Na^M = |('^)m| , the order of a modM , 

the subscripts a and M are suppressed except for in Sa and Sa,M- 
Sa=b = 1 if a = 6 , 6a=b = if a 7^ 6 . 

Definition 3.1. Let M G N, a G Z]^j. Then m = rria = mM = fTLa,M is the smallest 
positive integer with (A^^ = K*^)™!) 

m contains every prime factor of M, (18) 

4|mif8|M, (19) 

m = gcd(a^"' - 1,m) . (20) 

Moreover 

mo = gcd (m, M/m) , Mq = M/mrriQ = M/ gcd (m'^, m) . (21) 

In practice m is easily determined by Eq. (|28|). We begin with a technical proposition. 
Proposition 3.2. For Mi|M, a G Z]^^ we define 

f 1 if Ml odd, 

c= <^ (l + f + f) if Ml even, M = 4mod8 and m = 2mod4, (22) 

[ (l + f ) else. 

With this definition we obtain for = a'' 

c e Ki/n. , (23) 
a^™*^i^' = 1 + cmMiA;mod gcd (m=^Mi,M) , (24) 
3cfc G Z]^/^ with a^^'' = 1 + Cumk mod M , (25) 

E '^Mso+..a^c=omod^exp(^— Sofc + — s.a'^j, ifMlm^Mi. (26) 
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^ r /27ri 27Ti ,\ ^ 

^— Z^^mso+s,aHa^™-i)=omodMexp ( — Sofc+ — s^a 1, ifM|m. (27) 

Proof. Eq. (^31) is an immediate consequence of Def. 3.1. Moreover 

= {a"- - l) M,k (l + ^^4^ - l)) + (a--. - if'f ["f) {a"- - l)-' . 

First we show that the sum over j is a multiple of Mi. Obviously N 3 = 
(* • Since j has at most j — 2 prime factors for j > 3 and a^"" — 1 has by definition 

every prime factor of Mi we obtain N 3 (^^j^^^) (a^'" - iy~^k/j = {^^'Y) (^^"^ " lY'^Mi. 

The first term equals (a^™ - l)MiA; modm^Mi if Mi is odd, and (a^™ - l)MiA;(l + 
m/2) mod Ml if Mi is even. Finally the term M/4 can be added without harm if Mi 
is even and M = 4 mod 8. This gives Eq. (12^) . 

Let M2 = maXjgcd(A;,MJ). Then k2 = fc/Ma G Z^^ and with Mi = gcd(M2,M) 
we have gcd(m2Mi,M) = gcd(m2M2,M) = ^im^Ma + £2M for some £1,^2 e Z. With 
k^^ G Z^^, £3 G Z Eq. gives a^-'^ = 1 + cmk + iaiim^ M2k2k2^ mod M . We find 
Ck = c + i^iimk^' G Z^^„ from Eqs. (0), (||). 

Eq. (|26| ) is obtained if one plugs Eq. (|2^) into the definition ( |T0|) of g. We find 

(Ei 1 ^ /27ri , 27ri ^ 



/27r2 , 

= > exp Snfc H Snd" 



E E exp ( (A; + iV^MiA;') + ^s^a' (1 + cmMi^'; 



The sum over k' gives the Kronecker S and we get Eq. (]26|). 

Finally Eq. (^7]) will be obtained for Mi = 1 with Eq. (|36|) . This equation will next 
be derived independently of Eq. (pT]). □ 
Lemma 3.3. 



If m' is the smallest number that meets Eqs. (|IBp, (|TU]) then 

iV^ = iV^, , m = gcd (a^'"' - 1, m) . (28) 

, = A^2^opi...p, = Icm {N^ko , iVpi, . . ., iVpJ , 

^ Pi ""Pi 

if pi, . . are odd, prime and ko = min {k, 2) . (29) 

Nm = 1 m = gcd (a — 1, M) <;=^ a — 1 meets conditions 1 and 2 of Def. 3.1. (30) 

M = P"^, 2^ P prime: a is primitive, (a)pd = Z^^^ ^ m = P A = P - 1 . (31) 

M = 2'^, (i>3: a = lmod4 ^ A/'„ = 1 , a = 3mod4 ^ A^;^ = 2 . (32) 

M'|M^ Ar^^^,|Ar„ , mM'l"^- (33) 

rUm = m . (34) 
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{a)M = {a)m + mZM/m ■ (35) 

Mm ^ ' 

rrimo = mM/m ■ (37) 

^ = Mo, JL = !^_ (38) 

Nrno = gcd {Nm, NM/m) ■ (39) 

AT gcd (iV, A) m 

^--=gcd(iV.,A)' gcd(iV.,A) • 

Proof. Obviously m'|m and thus Njn'\Nm since (a)^' is a subgroup of (a)m- Since m is 

minimal we find Nm = Nm' and thus Eq. (^). Eqs. (^)- (l3^) follow directly from Eq. 



We use Eq. ( p6D with sq = and n = 1 to Fourier transform (a)^- We choose Mi 
with Sim? Ml = OmodM and obtain that g vanishes unless ca^simMi = OmodM. Since 
a G Z]^, c G 1^M/m (-^l- (HD) this implies SimMi = OmodM. Thus for non- vanishing 
g and M( = Mi/ gcd(Mi,m) we have sim'^M[ = OmodM which enables us to use Eq. 
(^) again (with Mi replaced by M() yielding the stronger condition simM[ = OmodM. 
Since m contains every prime factor of M and Mi|M we can continue until M{ - ' = 1 and 
finally we obtain 



g[N] {0,s) = VN/Nm6s, =0 modM/m9 

[iV,n](0,s). (42) 

This fixes (a)M to be {a^ + k'mmodM, k G IjN^^k' G Zm/to} (Eq. (p!3D). This is the 
statement of Eq. (^) and Eq. (|36|) follows immediately. 

From mM/m\M/m and mM/m\f^ (Eq- (0)) we get mM/m\f^Q- Eqs. (0) and ( ^3] ) give 
"^M/m = "^mM/„|"^mol"^M/»n, thus Eq. (|^ is Verified. 

Eq. (|3|) follows from Eqs. (|§), (§3): 

(§) MiV^„/,„ (§) MNm^^ (g) MiV^„ 0) iViV^„ 

JVAf/m — — — — • 

m mM/m rnmmo mm^ Nmirio 

To prove Eq. (R3) we begin with 



H ^ / NmrriM/m M 

^ gcd (AT^, A^M/™ j ^^^^ gcd (iV„, NM/m )^ gcd ( — 



'mo ^^rriM/m \ ^^"^Ml^ 

We will show that the right hand side equals m^. From Eq. (|25|) we know that there exists 

a c' G '^^M/,n with 



JU Nm ]\f ]\/f 

cm = a " — 1 = a — 1 = cmM/m^j mod — 
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Since c G '^M/m have also c, c" = cc' E TL^ui-m and the above equation imphes that 



there exists a k E 7^ with NmrnM/m/ Nmj^i^^ = kM/m + c"m. Finally we find for the right 
hand side of the first equation (c" G Eq. (|37|)) 

'"0 

gcd (mM/mA^m/A^m„/„, M/m) = gcd {kM/m + c"m, M/m) = nio . 

Since m\maA\a^'^"'aA — 1, we get NmlAN^ ^. On the other hand ^ smallest 
number with this property, so ^ = N^/ gcd(iVm, A). Moreover we find 



■m„A 



tea) 



AT^, .M 



N^M 



NrrM gcd (AT, A) (§) gcd (AT, A) m 



N,A gcd (N^, A) N,A gcdiN^,A)N gcd(Ar„,A) " 

Finally the first statement of Eq. (^) follows from Eq. ( ^0|) and Eq. (^). Moreover 

(g) gcd (at, A^M/^n) m Q A^M/^m (§)JI]) iW; 

gcd (^Nm,NM/m, 



m iv, 



mo 



mo 



□ 



Theorem 3.4. For M G N let = a'' modM. If a, s ■ Xq = s„ = J2]=i Sja^'^ G we 
obtain for \g\'^ defined in Eq. ( pil| ) and m, mo given by Def. 3.1, 

\g\'[N,X.]{so,s) = ^ E |^niV/iVMM,Xfc+Ar^,/^.](so,s-Xo) (43) 



mo ,.,2 



^lz^.„/iv.„, [^/^Af/m, ^fe+JV„,„.] (so, S ■ X, 



AT l-^" I^JVm/iVmn 

' mo " 



if 3A; G Z7v_ : So + s - (X7v^+fc-Xfc)/m = ^S^,^ mod mo , (44) 
else . 

Corollary 3.5. For all and all generators of the type 

Xk = cio!^ + C2k + C3 mod M (45) 

with period A^x we have 

|^nArx,M,X](so,s) = ^|^nAr',M',X](^^.o,^sj , with (46) 
M' = M/ gcd (cis„ M) , AT' = A^A,, , 

(H) apply. With m' = m^/, mj, = mo a/', 



and for the right hand side of Eq. (|46|) Eqs. 
Mq = M'/m'mQ we obtain 



|2 



AT; 



N^,N^,^M'^ '"'-w/'-™^ Mo 



[— —— , Xk+N,M'^,] (so, s • X, 



|^^A^,,X.] (so,s) = 



if 3A;gZ 
else . 



'■0 Nx 



+ 



M 



(47) 
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The denominators in the above condition are understood according to 

a/b = c/c/ mod I/Mq arfMg = cbM^ modbd . (4J 
Proof of the theorem. We start with Eq. ( pH] ) 



1 ^ /2Tci 2Tci fc / Afc , 

E e^p[-j^SoAk + —sy{a^'-l 



^ E E E exp('^soAk + ^Sa{a' + mk'){a^'-l] 

fceZjv„ fc'eZM/™ AfcGZjv ^ 

]^ E E '^s.(aAfc-i)=omodMMexp(^^soAA; + — s^a (a - Ij j 



^.a^ M/m ^ M " V 

[replace A/c by Ak — Ak' and A; by + NM/mAk'] 

1^ E E exp (^iV,,/„.o (AA; - Afc') + ^^aa'^ (a - a 



average over A/c' with Nm/tyi/N^^^^, 



which proves Eq. (^). 

Now we spht the /c-sum according to k = k^ + gcd{Njn, NM/m)k' = ko + Nm^k' . One 
can replace k' by NM/m/Nmok' since fc' G Z7v^/iv,„^ and NM/m/Nmo e Z]^^/^^^^. So = 

ko + NM/mk' and the A;'- dependence drops, since it amounts to a shift k k + k' of 
X-^ = a^+^M/mk (^cf. Eq. (|12D). The sum over k' gives simply Nm/N^o- 

Since mj^j^/^ = Omodm and m^Nj^j^^ = OmodM/m we have (?7i^iVM/m)^ = OmodM 
and Eq. (^) applies with a replaced by a^*^/™ and Sa by a'^Sa- With Eq. ( PD we get 



7710 



AT ^ N /N N 

//t0 



E ^^+,^„^+-M/™^V"-/-"-^"'"» -1)=0 mod M [-J^^M/mS^k + _ ,„a'=+^-/^'=' 



jv;; 



The /c'- dependence drops modM (Eq. (PBD) in the argument of the Kronecker 5. The 
fc'-sum interchanges with the Kronecker 5 and we obtain with Eq. ( pSD 



E '^^so+^iaa*(a^™"o-l)=OmodMl5'lz [^/^M/m,-^A;+Arjv^/^.]- 
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We use Eq. (^) with Mi = Mq, k = 1 yielding modM 



a 



NmMo 



M M /a^' 
1 = c— = — — 
mo mo V 



m 



m 



1 + -77^2\^ 



(i) M 



1 m 



mo V m 



Finally we divide the argument of the Kronecker b by M/rriQ and notice that, if 2|^^^, 
we have m/2 = mo/2 = — mo /2 mod mo- Since the /c-sum has at most one non- vanishing 
term we obtain the desired result. □ 
Proof of Cor. 3.5. With the definition s' = M'cis/M we find < = E"=i SjCi^'"^ = 
CiSa/ gcd{ciSa,M) G Z]^^,. With Nx = \cm{Nx,M) we obtain 



\9\ {so,s) 



1 
1 



^ I 2TTi 2Tii 
2^ ^^pIa^^o^ 



M 



2711 / Nx 



2 



27rz 



We can replace Ekez^^ by Nx/NxEk&j^^- This sum splits into Efcoez^/ Efe'ez^^^^, with 
k = ko + N'k'. The exponential depends linearly on k', thus the sum over k' gives a 
Kronecker 5, 



1 

N^ 



E ••• 



k&N^ 



N: 



X 



1 



E ••• 



The Kronecker 5 assures that s'q = N'{so/Nx + C2J2j=iSj/M) is an integer. The sum 
over k gives N'\g\'^[N' , M' , X^. = a^](so,s') and Theorem 3.4 applies. With integer and 
Eq. (EHI) the condition in Eq. becomes 



, f So C2Y.]=lSj\ M'ciSa fc 



A^' I — + 
\Nx 



M 



+ 



Mm 



ra' (a^™'-l) = ^52iM^modm'c 



-(52iM^ mod 



The expression in the brackets equals s ■ (Kx^,+k ~ X^). So we have confirmed the 
invariance of the condition in Eq. (^) as claimed in Eq. (^). Similarly 

\gnN'/NM>/m', M, X',_,x^^,^^^„] (4, s' ■ X'o) 

2 

.k+N,,„,,k' 



exp 27ri 

fe' 



-so + 



= l^nAT'/ATMyw, M', (AT'/AT;, ■ So, M'/M ■ s ■ Xq) 

= I^HA^x/iVA/'/m' , M, Xfc+A5,,^,/„^,.] (so, s ■ Xo) . 

Altogether we have proven that |5fp[A^x, M] is given by Eqs. (^31) , (44) and Eq. (0). With 
Eqs. dgg), (|3|), (H), (H) we find Eq. (g^). □ 
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4 Generators 

4.1 Xq = c, Xfc+i = aX/^modM 

We assume that 1 7^ a, c G Generators of this type are called multiplicative. The 
recursion formula implies Xk = ca^ mod M. The harmonic analysis of many multiplicative 
generators is easily obtained by Eqs. p^ ) and (^T]). 
Examples 4.1. 

1. Nm = l (Fig. 1). We find m = gcd(a -1,M), N = M/m, mo = gcd(m, M/m) = 
gcd((a - 1)2, M)/ gcd(a - 1, M) and Mq = M/ gcd((a - 1)^, M). 

With M' = M/gcd{sa,M) = M/sa,M we get = A^^, = 1, m' = gcd(a - 
1,M') and M/M^ = gcd(m'2 gcd(sa, M), M) = gcd( gcd((a - 1)^, M'^)^^, M) = 
gcd{sam'^,M) = mmogcd{sa,Mo) = mmoSa,Mo- Thus 

= Mo/sa,Mo if A^m = 1 . (49) 

We obtain from Eq. (|47|) 

I^P (So, S) = mogQ,Mo'^5o+es„(a-l)/m=lmos„,Mo<5 Mq modmos,,Mo ' 



In the case of a power of two modulus M = 2'^ with m = 4 (m = 2 is not possible 
for d > 2) we obtain for a proper choice of a (cf. Sec. |]) 

iVx = M/4, Qi = V2/4, Q„>2^M^/"-^ (51) 

If we look at correlations of k with Xk, disregarding higher n-tupels, we get A^^i = 
Ml = NQi = \/2M/16 reasonable random numbers with log(v^M/16) digits. 

A detailed discussion is postponed to the next section where we analyze the mixed 
multiplicative generator which is very similar but more popular than the multiplica- 
tive generator. 

2. M = P'^, P is odd, prime and a is a primitive element of l^pa, {a)pd = Tjpd, c G Zp^ 
(Figs. 2, 3). We find A/"^ = P - 1, m = P, = (P - 1)P^-^ mo = 1 if = 1 and 
mo = P iid>2. 



For the case d = 1 Eq. (|4^ ) is useless. However, we saw in Sec. ^]2| how to calculate 
\g\^ from the very definition. In fact it is even simpler to observe that is constant 
for all Sa ^ OmodP (Eq. (|12D). With Eqs. (|l|), (|15]) we obtain {P - l)6so=o + (P - 
l)\g\'^{so, Sa 7^ 0) = P with the result given in Eq. (|52D . 

If d > 2 we obtain for Sa ^ OmodP from Eq. (H \g\'^{so, s) = P/(P - 1) ■ Ek&p.^ 
Sso+saa''{aP-^-i)/PmodP = ^/(^ - 1) " ^so^omodP- For general Sa we use Eq. (||) and 
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obtain the result 



\9? (so,s) 


so = 


So 7^0 


Sa = 


pd-1 (p _ 1) 





X = d-1 


pd-1/ (p _ 1) 


p'5^=d-l/(p-l) 


X<d-2 





p'+\=,/ (P - 1) 



(52) 



where Sa,pd = gcd (s^, P'^) = , gcd (sq, P'^) = P^ 



We find that Qi = Qi(P^, P^) = v^(P- 1)/P is independent of a and even greater 
than 1. For n >2 only the case (sq, Sa) = (0, 0) contributes to Qn and the discussion 
is equal to the case with power of two modulus presented in Sec. |[ We find 

Nx = (P - I) P''-\ Q, = V2iP-l)/P, g„>2^P'/"-^ (53) 

From a mathematical point of view odd prime number moduli give good random 
number generators. In particular A^^i = P — 1, Mi = P whereas for a power of 
two modulus M we had = = V2M/16. So we need M > 16P'^/V2 to 
obtain power of two multiplicative generators which behave better than generators 
with powers of odd prime numbers. However, one has to take into account that 
computers calculate automatically modulo powers of two. Moreover, the power of 
two generator will be improved in the next sections until we achieve Qi = 1 (cf. Sec. 
|4^ ). As a byproduct a better behavior of Qn for n > 2 is obtained, too. 

Best performance allow prime numbers of the form P = 2^^' ± 1 0. In this case 
a ■ b = c\l^ + C2 leads to a ■ 6 = C2 =F cimodP. The extra effort, compared with 
a calculation mod2'^ is one addition and, which is more important, the calculation 
of c\. In Ex. 5.1 1 we discuss the improved generator with M = 2^^^ which can 
most easily be changed to M = 2^^^. Alternatively one may construct a generator 
mod 2^^^ — 1 which is a prime number. This generator will however be more time 
consuming and moreover it has worse quality A^i ~ 2^^'^, A^2 ~ 2^^'^, A^3 ^ 2^^'^, 
etc. vs. A'^i = 2^^^, N2 ~ 2®^'^, ^ 2®^, etc. for the generator with power of two 
modulus (cf. Sec. Sec. |^). So, power of two generators are more efficient than 



multiplicative generators with prime number modulus. The situation is slightly 
different if one considers multiply recursive generators with prime number modulus, 
analyzed in Sec. |4.5| and Ex. 5.1 2. 



Multiplicative generators with prime number modulus and primitive a produce every 
random number 7^ mod P exactly once in a period. We recommend to use a prime 
number modulus only if one needs this quality. 

3. M = 2^^, c = 1, a = 3 mod 8 (the cases a = 1, a = 5 mod 8 were discussed in Ex. 1). 
This example and the following ones are less interesting from the random number 
point of view. They are discussed to demonstrate how Eq. ( PI ) applies in less trivial 
cases. We restrict ourselves to Sa € and find 

, 8 if > 3 f 2^-^ if > 3 

2 if = 1, 2 ' I 2'^-! if = 1, 2 ' 
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mo 



8 if d > 6 

2 iid = 2 
1 if d = 1 



2'^'^ AT, 



mo 



2 if > 5 
1 if < 4 



The most complicated case is = 4 where we get |^p(so, s) = 6, 



so+Sa(a2-l)/8=0 mod2 " 



exp(^so) + exp(^(4so + asj) 
obtain the following table 



25,„=i^od2(l + cos(f (4so + (a - l)s,))). We 



M 


= 2'^ 




d 


> 7 


4 J2k=0,l ^so+Sa3'={a2-l)/8= 


=4 mod 8 


d 


= 6 


4Z]fc=0,l ^so+Sa3'={a2-l)/8= 


=0 mod 8 


d 


= 5 




26sQ=l mod 2 

Sa = 1, 3 mod 8 


Sa = 5, 7 mod 8 






So = 0,2mod4 








d 


= 4 


So = 1 mod 4 


2tV2 


2±V2 






So = 3 mod 4 


2±V2 


2tV2 






upper sign: a = 


-- 3 mod 16, lower sij 


j;n: a = 11 mod 16 


d 


= 3 




1 




d 


= 2 




2Sso=l mod 2 




d = 


= 0,1 




1 





(54) 



and Qi = V2/A, Qn>2 ~ M^/""^ for > 6 like in the case a = 5 mod 8. 
4. M = 10^, c = 1, a = 3, gcd(sa, 10) = 1. We find 



M = 


Nm Tu N mo 




\9\Hso,s) 


d>9 


4 80 5 ■ 10'^-2 80 


4 


20 l]fc=0 ^so+3*Sa=40 mod 80 


5<d<8 


4 80 5 ■ lO-^-^ 5 . 2^- 


-4 4 


5 • 2*^ ^ X]fe=0 '^so+3'=Sa=0 mod 5-2^-* 


2 < d < 4 


4 5-2'^ 4-5'^-i 5 


4 


5/4 ■ 5so^0mod5 


d=l 


4 10 4 1 


1 


1/4 + ^so^o mod 4 








(55) 


and Qi = V2/20, g„ < g„(0, -3, 1,0,.. 


,0) = 2- 


10=^/2-d^ if ci > 8. 


M = 2^ ■ 7, 


c = 1, a = 3, gcd(sa, 14) = 


1. This example illustrates all aspects of 



calculating We find A^^ = 6, m 



■ 7, 



3, M/m = 2^ 



mo = 8, A^„ 



M/m 



l)/m = -3 mod mo. Eq. (|4|) gives |^p(so,s) = 4/3 



Efe=o,i^so-3^+^s.=4mod8l^ll3[A^ = 24,M = 27-7,a = 34](so,3'=sJ. 
Firstly, = 3'=(1 - 4)s„ = 3^Sa + 4 mod 8. 

Secondly, in the we sum over k' G Z3. Since 8 G Z3 we can replace k' by 8k' 
yielding |^ | ^ [AT = 3, M = 2^ ■ 7, a = 3^2] (^o, 3'^s J . 
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Finally, we notice that 3^^ = lmod2^. With the Euclidean algorithm we find 
integers ci, C2 with 1 = 2^ci + 7c2- This plugged into the definition of g gives 
exp(^Sa332fe'(^2^ci + 7c2)) = exp^^Sa^^^'"' ci + ^SaC2). The second term is inde- 
pendent of k' and drops in Since 3'^^ = 2 mod 7 and ci = 2~^ = 4mod7 we 
obtain with Eq. (|T2|) 

|^|^(so,s) = | 5^„+3fc^^=omod8l^n^ = 3,M = 7,a = 2] (so,3''Sa) . 

fc=0,l 

An explicit calculation (cf. Eq. (PT])) is needed to obtain 



I^P (^5o,s) 


Sa = 1,9, 11 mod 14 


Sa = 3, 5, 13 mod 14 


So = mod 3 


8/3 


8/3 


So = 1 mod 3 


2/3 ■ (7 - V21) 


2/3 • (7 + V21) 


So = 2 mod 3 


2/3 ■ (7 + V2I) 


2/3 ■ (7 - VWj 


|^P(so,s) = 


= 





if 3A; e 0, 1 with 
So + 3''Sa = OmodS 



else . 



(56) 



4.2 Xo = 0, Xk+i = aXk + c mod M 

Let 1 7^ a, c G Z^.^. The recursion formula implies Xk = c{a^ — l)/(a — l)modM. 
Generators of this type are called mixed multiplicative. 

There are two ways of deriving the mixed multiplicative generator from the multi- 
plicative generator. 

Firstly, we expect from a good generator to have random differences AX^ = X^+i— X^. 
For the multiplicative generator we have AX^ = (a — l)XfcmodM and in deed for the 
case M = P'^, (Ex. 4.1 2) where the multiplicative generator was good, the difference 
generator has the same quality as the original one. However if M = 2*^ obviously a — 1 
and M have common divisors. If e.g. a = 5 mod 8 then {AXk)k produces only multiples 
of 4 which deteriorates the randomness of the sequence. This suggests to use the sum 
generator SX^ = I]j=o -^j = '^(a'^ — l)/(a — 1) which is the mixed multiplicative generator. 

Secondly, we observe that a — 1 divides — c = c{a^ — 1) modM. Since Xk — c and 
Xk differ only by a shift they have the same quality. The flaw that Xk — c produces only 
multiples of gcd(a — 1, M) can be corrected by a division by a — 1 which yields the mixed 
multiplicative generator. 
Theorem 4.2. Let B be deflned by 

B = max fc gcd (a - 1 , M'^) (57) 

then the Fourier transform of the mixed multiplicative generator Xq = 0, Xk+i = aXk + 
cmodM, 1 7^ a, c G l^^j is given by Thm. 3.4 and Cor. 3.5 with MB and XB instead of 
M and X, respectively. 

Proof. From a,c, {a — 1)/B G Z'^^ we conclude a,c, {a — 1)/B G Z'^.^^- If we denote the 
inverse of (a - 1) /B mod MB by ^ we have Xk/M = c-£^{a^ - I) /MB mod 1. Cor. 3.5 
applies with XB and MB instead of X and M. □ 
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Example 4.3. 

1. NruMB = 1 which is equivalent to 

1. a — 1 contains every prime factor of M, 2. 4| (a — 1) if 4|M. (58) 
(Figs. 4-8). With 

6= gcd(a- 1,M) = gcd(5,M) (59) 

we find niMB = gcd(a - 1,MB) = B, ttiomb = b, Mqmb = M/b, Nx 
MB/niMB = M and s • (Xi - Xq) = s^c. We obtain from Eqs. (0), (||) 

I^P (so, S) = bSa^M/b5so+csa = hbs„„M;HS m mod • (60) 

'"'a,M/b 

In particular for sq = we observe that g is only non-zero if M/bsa,M/b is odd and 
Sa = Om.odbsa,M/b- From the second identity we conclude Sa = OmodM and thus 

|^P(0,S) = M6s,=0modM (61) 

For n = 1 we get Sa = Si and g is the Fourier transform of the uniform distribution. 
So, Xk has the property to produce every number in Zm exactly once in a period. 
The generator is called to have a full period. 

We proceed with a discussion of the parameters a and c. 

(a) Choice of c. We can restrict ourselves to 1 < c < 6/2 since every c emerges 
from a c G (0, 6/2) via translations (X^ X^+Ak — ^Ak = ci^^^k) or reflection 
(Xfc 1-^ —Xk). We can determine c by the condition that there should be no 
small (so,Sa), gcd(M, Sa) = 1 with sq + SaC = 6/2 ■ 62\M/b^odb. If M/b 
is odd c ~ Vb gives Qi(so ~ \/b,—l) ~ -\/6 + 1/6 ~ 6^-*^/^. In the case 
where M/b is even and 6 is small the choice c = 1 is best with the result 
gi(so ^ si ^ 6/4) f« ^2 ■ (6/4)/6 = ^2/4 ^ 0.35. 

In the case of Fig. 5 with the 'wrong' choice c = 3 one has a = 9 mod 16, 6 = 8 
and the smallest (sq, Si) with non-vanishing g is (1,1). Since |^P(1, 1) = 8 we 
obtain Qi{l, 1) = V2/8 ~ 0.18 (notice the correlations perpendicular to the 
(l,l)-direction in Fig. 51). With the right choice c = 1 (Fig. 6) it takes an 
(so,si) = (1,3) (or (3,1)) to get = 8. Therefore Qi(l,3) = VTO/8 ^ 0.40 
which means that the random numbers are more uniformly distributed in Fig. 
61. The large value of Qi{l, 3) is yet misleading since Qi = Qi(— M/8, M/8) = 
V2/8. However Ig]"^ assumes the small value of Qi at much less sites as in the 
case of c = 3 which means that the choice c = 1 is better than c = 3. 

Notice the similar pair distributions in Fig. 511 and Fig. 611. In general the 
quality dependence on c can not be obtained by the standard spectral test 
(corresponding to sq = 0) since the n-tuple distributions are only shifted by 
changing c. 
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(b) Choice of b. In general b should be as small as possible in order to prevent 
I^P from being concentrated on too few points. If M/b is odd, c ~ \/b then 
Qi behaves like If M/b is even, c = 1 then Qi{so,Si) = Qi{sq ^ si ^ 
6/4) « V2/A for small (sq, si) (cf. (a)). However Qi{-M/b,M/b) = V^/b 
which forbids large values of b. 

For a power of two modulus the smallest value possible is 6 = 4 which implies 
Qi = -\/2/4 (Fig. 4). In particular M = (Fig. 7) should be avoided since 
in this case b > 20. 

These arguments require sq ^ 0. They are not obtained by the standard 
spectral test. 

(c) Choice of a. Up to now we have evaluated 15^^(50, Si) for n = 1 which is given 
by b and c. In order to determine a more precisely we have to look at the 
distribution of n-tupels for n > 2. In this case Q„ = mmsQn{so = 0,Sa = 
0) = mins:s„=o|s|/M. A further discussion of the choice of a is postponed to 
Sec. ^. We will see that a reasonable a gives Qn ~ M^/"~^. 

Since Sq = this part of the choice of a is identical with the standard spectral 
test. 



Let us summarize the result for M = 2^ 



c=l, a = 5mod8, maXamins:s„=o|s|, for n = 2, 3, . . . gives (62) 
Nx = M, Qi = V2/A, Q„^Mi/"-^ (63) 

A loss of randomness for n-tupels is avoided if one takes gcd(ri, M) = 1. 

2. a = 3 mod 8, M = 2'^. If Sa,M < M/32 then B = 2, tji'mb = = 8, ^o'ms = 

M/32s,,M, Nm'^^^ = Ar„„,^^ = 2 and from Eq. (g^) we get 

1-12/' ^ _ 8NxSaM r /^/|\ 

\9\ [S0,S)- — 2^ M ^^_^_3fc^a+i^^^4^^^^^^ ^^^jg^^_^^ . (^D4j 

^^-^ fc = 0,l ^ ' »a,M 

We obtain A^^ = M/2 (cf. Eq. (|l3D) and either gi(l,l) = ^2/4 or Qi(-l,l) = 
v/2/4. The period is doubled compared with the multiplicative case and the quality 
remains unchanged. Since a = 3 mod 8 has only half the period of a = 5 mod 8 the 
latter should be preferred. 

The only improvement we achieved was an increased period by a factor of B. The 
mixed multiplicative generator mod M equals the multiplicative generator mod Mi?. For 
a power of two modulus we still have Qi = a/2/4 < 1. 

This does not mean that one can not use the mixed multiplicative generator if the 
modulus is large enough. In Fig. 8 we see that all problems dissolve if the modulus is 
much larger than the range of random numbers used. Nevertheless it will be profitable to 
look for improvements in the next sections. 

We close this section with a well known theorem which is a corollary of our general theory. 
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Corollary 4.4 (Greenberger Hull and Dobell p]). The mixed multiplicative generator 
has a full period if and only if (|58[) holds. 

Proof. We saw in Ex. 4.3 1 that (|58D implies a full period. If on the other hand NmMB > ^i 
due to Eq. (p9D there either exists a prime factor P ^ 2 of M with A'^p > 1 or 4|M and 
A^4 > 1. In the first case we look at the generator modP and find B = 1. Hence 
Nx mod p = Np < P and the generator has no full period. In the second case a = 3 mod 4 
and Vfc: Xk = 0mod4 or Xk = cmod4. □ 



4.3 Xq = 0, Xk+i = aXk + ck mod M 

Let 1 7^ a, c G "Z^j. The recursion formula gives X^ = ^^(^r^^ — k) modM. 

The generator of this section improves the mixed multiplicative generator according 



'mm 



to both ways presented in the last section. Firstly, it is the sum generator J2jZo -^j 
X)j=o c(a'' — l)/(a — 1) of the mixed multiplicative generator. Secondly, we observe that 
X™™ — ck = cX]j=d(a-' — l)modM can be divided by a — 1. On the other hand the 
randomness of X^ should not differ much from that of Xk — ck and the division by a — 1 
improves the generator. 

Theorem 4.5. Let B be defined by Eq. (^) then the Fourier transform of the generator 
Xq = 0, Xk+i = aXk + cfcmodM, 1 7^ a, c G Z]J^^ is given by Cor. 3.5 with 

X ^ XB^ , M ^ MB""' . Nx = 1cm (^Nm, 2^^\^' M) (65) 

is the period of X. 

Proof. Analogously to the proof of Thm. 4.2 we find that X^/M = c((^)^(a^ ~ ^) — 
■^Bk) /MB^ mod 1 where ^ is the inverse of (a - l)/5 in Zmb^ ■ With Cor. 3.5 it only 
remains to show that the period of X is given by Eq. (|65|) . 

Since = Xi = Xjvy+i = aX]^^ + cNx = cNx modM we know that M\Nx and since 
d^x _ ]^ _ OmodM implies Nm\Nx we have 1cm (X^,, M)\Nx. If we split M into its prime 
components M = 2'^p\^ ■ ■ ■ pl\ it suffices to show that X^^ = Omodp^'' Vj = !,...,£ 
and X2k = 2'^^^mod2^. If Np. > 1 then a — 1 G T^p. and one obtains immediately 

Xnx = ^(^^^^ - Nx) = OmodpJ^ The case Np^ = 1 is treated in Ex. 4.6 1. For 
M = 2 we get X2 = 1 and for M = 2^ the result follows by induction since 

X2M = (^—^ - 2m] = ^— ( - M 1 + c ( ) = 2XMmod2M 

a — lya — 1 J a — lya— 1 

for M > 2. □ 
Examples 4.6. 

1. Nm^^g2 — 1 which means that a — 1 contains every prime factor of M and 4|(a — 1) 
if M is even (Fig. 9). We find niMB'^ = ^omb'^ = B, Mqmb'^ = M and for odd M 
we obtain 




I - 12 / \ _ Nx (- x^mm _ ^ ^ 

\9\ 1'50, Sj - — S<j,M(J^so+s-X™=Omods,,M ) " C ^ _ ^ • (Obj 



4 GENERATORS 



21 



With Eq. (|T^) we can conclude that Nx = M since any Nx = AM leads to 
|^P(so, s) = unless A|so. This completes the proof of Theorem 4.5. 

If M is even we have Nx = 2M yielding 

I^P ("So, S) = 2Sa,M5so+2s-Xff'"=s,,A/5^^u_ mod 2s„,m • (67) 

For n = 1 we have Sa,M = si^m = gcd(si, M) and s • X™™ = siX™™ = 0. Thus for 
odd M we have Qi = Qi(0, si,a/) = 1- If M is even then Qi(l, 1) = ^2/2 and if 
3fc e N so that 1 7^ M/2'= is odd then Qi(0, 2'') = 1/2. Therefore 

r 1 if M is odd 

Nx = 2^^f^m, Qi = \ V2/2 i{M = 2'^ , Qa ~ M-^/^ Q„>3 ^ M^/^^-^^-^ 

[ 1/2 else 

(68) 

In case of a power of two modulus we have increased Qi from \/2/4 for the mixed 
multiplicative generator to V2/2 ^ 0.70 and Nx from M to 2M. Notice also the 
better behavior of Qn for n > 2 which we will derive in Sec. ^. 

For So = we obtain \g\'^{0,si) = 2^^\^^ si^m^m/si^m odd- Only for M = 2'^ this is 
NxSsi=o, the Fourier transform of the uniform distribution. In this case one obtains 
two full periods, every random number occurs twice in a period. Analogously to the 
proof of Cor. 4.4 one can show that for N^.^^^^ > 1 a full period is not possible. 

Thus the generator Xk+i = aXk + ck mod M produces exactly uniformly distributed 
random numbers if and only if M is a power of two and a = lmod4. 

The distribution of n-tupels is no longer a lattice as can be seen from Fig. 911. So 
even for the standard spectral test (sq = 0) one would need Eq. ( ^T]) to analyze 
the generator. We skip a further discussion here since in the next section we will 
present an improved generator for a power of two modulus. 

2. M = 2^=^, a = 3 mod 8. Let Sa,M < M/IQ then B = 2, m^.^^a = m'oj,,^^2 = 8, 
^0MB2 = M/16sa,M, ^m'^jg2 = ^<mb2 = ^ and we obtain 

\g\^ (so, S) = 8S,,M E '^^o+s.(X— +X— )=8s.,m53,|^ modl6«.,M " ^^^^ 

k = 0,l ' ^a,M 

For e.g. c = 1 we obtain Qi(3, 1) = VTO/8 ^ 0.40 or Qi(l,3) = v^/8. The 
improvement in the quality is very little in comparison with the mixed multiplicative 
generator. The period increases from M/2 to 2M. The generator with a = 5 mod 8 
behaves clearly better. 

3. M = P"^, P is an odd prime number and a is primitive, {a)pd = Zp^j (Fig. 10). We 
find 5 = 1, Nx = (P — 1)P'^. In the following ^ is c times the inverse of 
a — 1 mod P"^. We have to distinguish three cases: 

1. Sa pd = OmodP*^, m' = tjiq = Mq = 1, Nm' = N^' = 1- The condition in Eq. 
(|7p is So + (P - l)s ■ (Xi - Xo) = Omod(P - 1)P'' ^ so = OmodP - 1 and 
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so + (P - 1)C{S, - - 1) = SO - (P - 1)^ E ■=! S, 

the pre-factor (P — 1)P'^ we obtain the first hne in Eq. ([70|). 



OmodP'^. With 



2. s^^pd = P'^-\ m' = P, m'o = 1 = M^, A^"^/ = P - 1, A^^j, = 1- We obtain from Eq. 



(|7|) |^|2(so,s) =P^|^||^_J(P-l)P^X.](so,s)5,„+,.(x,_,-Xo)=omodP^. Since P^G 
Z^_i we get P'^|^P[P-l,Xpd.](so,s)(5,^+_^(,^(^P-i_i)/(^_i)_(P_i)^n^^,.)^o^„dpd = 

P'^|^P[iV = P - 1,M = P,X, = (^a^'''^](so,s/P'^-i)5,„=(p_i)_l_^.^^,^^„,p.. 

Since a^"^^ = a^modP we can use the result of the muhiphcative generator, Eq. 
(|52D, d = 1, yielding the second line in Eq. (|70|). 

3. Sa,p<i = < P'^'^ m' = rn'o = P, = p^-^-s, A^^, = A^^, = P - 1. We 
obtain \g\\so,s) = P^+V(^- 1) -EfceZp., 5.o+s.(Xp_i+,-X,)=omodpA+2 = P^+^/{P- 

5,o__^(P_i) ^"^^ .,=fc's,p mod pA+2. Altogether, 



I^P(so,s) 


fi = d , So = OmodP — 1 


fi = d , So 7^ OmodP — 


1 fi < d 












A = 


P7(P-1) 


p'^+V(p-i) 





A < d-2 








P^+25,=A+l/(P-l) 



With s„,p. = gcd P'^) = P^ , gcd ((a - 1) So - c (P - 1) E]=i s„ P") = P^ 

(70) 

We find that a proper choice of a and c leads to Qi ~ \fP{P — 1)/P^ ^ P The 
by a factor of P longer period compared with the multiplicative case (Ex. 4.1 2) is 
to some extend compensated by a worse relative quality Qi. 

If we set So = and n = 1 we find A = /i and therefore \g\'^ = NxSst^=o- The 
random numbers are exactly uniformly distributed, the generator produces P — 1 
full periods. 

Although NxQn is larger than for the multiplicative generator we can not recom- 
mend the generator. If one is not interested in having full periods, generators with 
a power of two modulus are more efficient. If one needs a full period the velocity of 
the generator is not important since M has to be comparatively small. Then only 
the relative quality Qn is essential and the multiplicative generator is better. 

4. Xfc = c(fc^ — k)/2. This quadratic polynomial is obtained with a = M + 1. Since, 
of course, Nm^^^^ = 1 Eqs. (^) and (^) of Ex. 1 for odd and even M are valid, 
respectively. Although this generator has reasonable Qi it is too simple to have 
good pair and n-tupel distributions (cf. Sec. 



We close this section with a lemma which will play an essential role in Sec. ^A . 
Lemma 4.7. Let A e M' = MB'^ /sa,MB^, fn' = rriaA^M', i^'o = ^OaA,M'-> ^ 
M'/m'mQ. Moreover M' = MP^/s^jQ-p2, m' = tti^a m^, rh'o = fnoa^^M', = M'/rn'mQ. 
If one of the following conditions holds there exists at most one C G Tja with ^[A^, M, X, = 
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{XA,+C+j-l)j=l...n] 7^ 0. 



1. 3Ma e N with Nm^ = A and Ma\ {M' /m) . (71) 

2. 3M\M with a = ImodmQ and 

^^A.o + s . (X^.^, - Xo) = ^5,M^ mod ^ . (72) 

Proof. Since I (M'/m') we get Ma I M'. Thus Ma| gcd(a^-l, M')| gcd(a^^'"'-l, M') = 
m', hence M^lmg. Because of the condition in Eq. ( PI) and Eq. (|i6|) with Xi?^ and MB'^ 
instead of X and M we obtain as a necessary condition for ^[X, = {XA»+Ci+j-i)j\ 7^ 7^ 
g\X, = {XA,+c2+j-i)j]- 

3ki, k2 , with 5Z (-1)' II (^A(7v„,+fc,)+c,+i-i - XAfc,+c,+i-i) = Omodmo 

^ _sa cB' a^""-' - 1 fAkM _ ^Ak,+cA ^ Omodm[, ^ a^^ = a^^ modM^ 

and therefore Ci = C2. 

In case of Eq. ( |72D we may relax the condition in Eq. ( ^71) and replace (after multipli- 
cation with M) M/Mq by M/M^. Since a = ImodmQ we see from the above calculation 
that in this case the condition is independent of C. Therefore all the ^[X, = {XA»+c+j-i)j] 
are zero if ^[X, = {XA»+j-i)j] is zero. The dependence on ki, /c2 drops, too, and we end 
up with the condition in Eq. (^). □ 



If g is decomposed ^Phg[Nx] = Ecez^exp{2mCso/Nx)^JNx/Ag[Nx/A,{XA.+c+J-l)J] 
only one of the terms on the right hand side can be non-zero. Therefore XxIfi'Pi^Js:] = 
Ec&, Nx/AmNx/A, (Xa.+c+.-i).]. 

This shows explicitly that the quality of n-tupels deteriorates if A\n and A\Nx because 
the period reduces by A and maxceZA \g?Wx/A, {XA,+c+j-i)j] = A\g?[Nx,X,]. 



4.4 Xo = 0, Xfc+i = aXfc + cint(/c/t)modM 

Let 1 ^ a,c e Z^, 2\t\M and Nmb^ = 1 (Fig. 11). 

This section is restricted to even M since for odd M the random number generator of 
Sec. was already satisfactory. We only discuss Nmb^ = 1 since this suffices to find a 
power of two generator with Qi = 1. 

If one tries to improve the generator of the last section by building sum generators one 
would obtain the recursion relation Xk+i = aXk + P{k) modM, where P is a polynomial 
of degree greater than one. Such generators may be interesting from a mathematical point 
of view, however they are not easy to analyze and the calculation of the random numbers 
becomes more time consuming. 

We follow the second way to improve the generator of the last section. Looking at this 
generator mod 2 one obtains the sequence 0, 0, 1, 1, 0, 0, 1, 1. . . and to make all random 
numbers even one could perform the map X^. 1-^ Xj. + mt{k/2) which does not seriously 
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affect tlie quality of the generator. If we afterwards divide the random numbers by 2 
we should end up with an improved generator. A slightly more general framework leads 
to the recursion formula of the generator of this section (cf. Prop. 4.9). It is perhaps 
surprising that even this generator can be Fourier transformed algebraically at almost all 
points (so, s). 

It would be possible to analyze this generator for Nmb^ > 1- However this would 
cause some extra effort and for practical purposes this is not needed. 



The Fourier transform of the generator is again derived from Eq. (47). 
Theorem 4.8. Let Xq = 0, Xk+i = aXt + cint(fc/t) modM, I ^ a,c e Z^j, 2\t\M and 
= 1 then 

\g\^ (so, s) = 5'5so+s-X--=omod^„,M with S = Sa,M if Sa,M\M/t . (73) 

X™™ = c(a^ — l)/(a — 1) is the mixed multiplicative generator related to X. The period 
is Nx = Mt. 

Proposition 4.9. There exists a c' G l^^j, a G Z, odd, and a periodic function R, 
R{k) = R{k + 1) G Q, RiO) = 0, so that we get with X'^ = ^{s^ - k), 

Xk = ^(^Xl-^k + R{k)^ modM . (74) 

Proof. The equation is trivial for k = 0. We can assume by induction that the equation 
holds for k and prove it for k + 1. With k = tk' + k^, < k^ < t we get X^+i = 
^^{aXl-ad'k/2 + aR{k) + tcmt{k/t)) = \{X'i^^^ - c'k - ad'k/2 + aR{k) + c{k - ko)). If we 
compare this with i(X[+i -(i'(A; + l)/2 + i?(A; + 1)) modM we find (a) c' + ad'/2-c = d'/2. 
Further on, R{k + 1) is determined by R{k) and ko. Therefore the periodicity of R 
is equivalent to R{t) = R{0) = 0. This means (b) Xt = ^(^^^ - 1) - d'/2. To 
solve Eqs. (a) and (b) we need the following observation. Due to Eq. ( pSj ) (replace M 
by MBt"^, B = m = gcd(a — 1,MB)) we know that there exists a q G I^Mt^ with 
a* — 1 = CfBt mod M Bt"^ . Thus ^^riy ^ '^Mt has an inverse modMt which we will 

denote by ^^°tZi ■ Now we can solve Eqs. (a) and (b) for c' and d' modMt, 

, t(a-l),, ,,,, , 2c t(a-l)\ J 

a* — 1 a — ly a* — ly a* — 1 

Since every prime factor of M is contained in a — 1 and c G we have {a — l)Xt + c G Z]^ 
and therefore c' G Z]^. It remains to show that d' is odd. We obtain mod 2: 



,, a — 1 c I a — 1 \ a — 1 c^^lTl a — It 
d' = 2- t = 2- Xp'^- ^ = 2- = 1 mod2 



a* — la — l\^a — 1 J — 1 ^ a* — 12 

□ 

Proof of Thm. 4.8. We start with Prop. 4.9. Let k = tk' + ko and |P(A;o)P = 1 then \g\^ 
equals 

2 



1 



E P (ko) E exp (tk' + ko) + ^ E f^;.'+.o+.-i - f (^^' + ^o) 
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Now we combine the d'tk' /2-teTm with the sotk'-teim yielding (sq — Nx/2Mt ■ d' J2]=i Sj) ■ 
{tk' + fco). For the moment we set Nx = 2Mt since we want to apply Lemma 4.7 on the 
generator X'. Later we will show that actually we have Nx = Mt. 

We assume that one of the conditions (|7TD , ( [72D is valid (which we will show later). 
Then only one term in the sum over ko can be non-zero. Therefore \g\'^ = 1/Nx ■ 
X]fco l-P(^o) Z]fe'... P P{ko) drops. Afterwards we interchange the sum over ko and 
taking the square once more and recombine the fco-sum with the fc'-sum. We end up with 
the squared Fourier transform of X' which was already analyzed in Ex. 4.6 1. The result 
(Eq. (13)) is 



\9\^[Nx - 2Mt] (so, S) - 2Sa,Aft5so-rf'X;"=i«J+2s-X;™=s„,Mi<52|^wj_ mod2s 



We find -d' E"=i Sj + 2s ■ X^-- = -d' E"=i sj + 2c'{sa - E"=i Sj)/{a - 1) = (d'(a - 1) + 
2c') (sa — Z)j=i Sj)/{a — 1) — d'sa and with Eq. (a) in the proof of Prop. 4.9 this equals 
2s • X™™ — d'sa- Now we assume that Mt/sa,Mt is even. Then Sa/sa,Mt is odd and since 
d' is odd d'sa = Sa,Mt^od2sa,Mt- Altogether the equation in the Kronecker 6 in \g\'^ gives 
So + 2s • X™™ = d'saSMt/s^^Mt odd ^od2sa,Mt- We recognize that Mt/sa,Mt can only be 
odd if Sa is even so that |^p[iVx = 2Mt](so,s) vanishes for odd Sq. From Eq. (p^) for 
ci = C2 = 2 we see that X has actually the period Nx = Mt and since d' is odd 

\g\^{Nx = Mt] (so,s) = Sa,Aft5so+s-X--=f 5Mt/.„_^,, odd mods„,Mt • 

Now we check the conditions of Lemma 4.7. We set Mt = max^gN gcd(ti?, t^). Since 
rriMt = gcd(a - 1, Mt) = Mt/t there exists a G Zt (Eq. (|5D) with a'' - 1 = 



CkkMt/tmodMt, hence Nmi = t. Assume Sa,M\M/t. Since Mt consists only of prime 
factors of t and gcd(Mt,t^) = gcd(i?t, t^)| gcd(M_B/sa^M, ^^) it suffices to show that 
gcd{MB/sa,M,t'')\gcd{M'/m',t'') for all k (M' = MBH/sa,MB'H)- Since Sa,M\M/t we 
get gcd{sa,M,t'') = gcd{sa,MBH,t'') and therefore gcd{MB / Sa,M,t'') = gcd{M' / Bt,t^). 
Finally from Eq. (H) we get gcd(a* - 1, MBH) = Bt, hence m' = gcd(a* - 1, M')\Bt. So 
Eq. ([fl|) holds if Sa,M\M/t. Moreover, in this case Sa,Mt = Sa,M and 2\t\Mt/ Sa,M which 
proves Eq. (|73D if Sa,M|M/t. 

If Sa^M does not divide M/t we take M = tSa^M in Eq. ([721). So M\Mt (notice that 
M has to be replaced by Mt in Lemma 4.7). Since M' = Sa^uB'^t/ Sa,MB^t\BH and 
m' = gcd(a* — 1, M') = gcd(i?t, M') we have mo|(M'/m')|i?. Therefore a = Imodrfig 
and = M' /m'm'^ = 1. We can use Lemma 4.7 if ^t(so - d' E"=i Sj) + s ■ (XJ - X[,) ^ 
^<52|M^modM. Nows-(X[-X'o) = -^^(s^^l^ -t^U s,) = + t^(s, - E •=! «,) 

and with Eq. (a) in the proof of Prop. 4.9 we get ^(sa — I]"=i Sj) = s ■ X™™ — |-(sa — 
X]"=i Sj). Collecting all pieces we obtain the condition tso/2 + ts ■ X™™ + Sa{X^ — d't/2) ^ 
^^2|M^modM. 

Since rUat^MBH = ^Oa\MB^t = Bt we get with Eq. ( ^9]) Mq = M/ gcd{tsa, M). There- 
fore 2^ = tgcd(ts„,M)/2 = OmodM if 2|M^. With Eq. (0) we get s„(X; - d't/2) = 

SatXt = mod M. We can divide the condition by t and get so/2 + s ■ X™™ ^ mod Sa,M- 
This condition holds for odd sq which completes the proof that g[Nx = 2Mt]{so,s) 
vanishes for odd Sq. The period is A'x = Mt and due to Eq. (|T3p we may replace So/2 by 



4 GENERATORS 



26 



So if we calculate g[Nx = Mt]{so,s). This means that Eq. (|73D is valid if the Kronecker 
6 vanishes. Otherwise the equation holds trivially with S = \g\'^{so,s). □ 

If Sa,Ai fails to be a divisor of M/t then remains undetermined for Mt/sa,M values of 
sq. At these points we only know that the average of is Sa,M (Eq. ([T5|) ). 

Therefore one should choose t as small as possible, namely t = 2. In this case 
is known if Sa,M\M/2 and an explicit calculation using Eq. ( fflD (note that a'' — 1 = 
(a^ - l)k/2 - (a - odd /2 mod 2^2) gives 



S = M\^ + cos\^—\^So + 2cJ2 ^2_i ^jjjj ifs, = OmodM, t = 2. 

(75) 

In particular for n = 1 we find |5'p(so,0) = 2M6so=o (cf. Eq. (P^), and n = 2 gives 
|5p(so, Sa = 0) = M(l + cos(7rso/M)) ^ 2M for small sq. 

If moreover M is a power of two then is completely determined by Eqs. (0), ([75|) 
and the generator can also best be implemented (cf. Ex. 5.1 1). We obtain 

Nx = 2M , Qi = l, ~ ^ ^ j^^i/(n-i)-i for n > 3 (76) 

as will be shown in Sec. ^. It is advantageous to have two parameters a and c at hand 
to optimize the quality of higher n-tuples and not only a as in the case of the mixed 
multiplicative generator. 

The generator does not provide full periods since |5'p(0,si) = si^m 7^ 2M6si=q. How- 
ever the deviation from an exact uniform distribution is not larger than in a finite true 
random sequence. If one does not use the entire period of the generator (and this is not 
recommended because of the n-tupel distribution) the feature of having a full period is 
anyway irrelevant. If, for some reasons, one insists in a full period we recommend to use a 
multiplicative generator with prime number modulus or the multiply recursive generator 
which will be analyzed next. 

The generator of this section behaves in every aspect better than the widely used 
mixed multiplicative generator. This is also confirmed by the figures (cf. Fig. 4 vs. Fig. 
11). The extra effort in calculating random numbers is little (cf. Ex. 5.1 1). If ra-tuples 
are used one should take odd n and occasionally omit one random number. 

Nevertheless the most essential step for producing good random numbers is to use 
large moduli (cf. Fig. 8 and Sec. 



4.5 Xo = l, X_i = ... = X_r+i = 0, Xk+i = ar-iXk ^ ... ^ aoXk-r+imodP 

Let P be a prime number and Xk have the maximum period of P*" — 1 (Fig. 12). In 
this case the random number generator has a full period in the sense that every r-tuple 
(Xo, . . ., Xr_i) 7^ (0, . . ., 0) occurs exactly once in a period. 
Theorem 4.10 (Grube 0). Let 

P (A) = - a,_iA'-^ - . . . - ao. (77) 

The corresponding generator has maximum period if and only if P is a primitive polyno- 
mial over Zpr. 
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The proof is found in Satz 2.1]. 
Theorem 4.11. Let Xk be defined as above and = (s-Xfc)o<fc<r = (I]j'=i SjXk+j-i)o<k<r 
then \g\'^ is given by the following table. 



I^P(so,s) 


so = 


So 7^0 


Sa = 


pr _ I 





s„^0 


1/{P'- - 1) 


pry (pr _ ^) 



Proof. First we notice that (s ■ Xfc)^ obeys the same recursion relation as {Xk)k since 
s ■ Xfc+i = J2]=i SjXk+j = J2£=i Q-r-i Yj^=i SjXk+j-i = Y^i=i ttr-e s ■ Xfc+i-^. So, (s ■ Xfc)fc is 
either identically zero or it has maximum period. In the latter case every number G Zp 
is produced P^~^ times in a period and the zero is generated P^~i — 1 times. Since the 
same holds for (s • (X^+Afe — Xfc))^ we get 

/27ri \ 

( ■ (Xfe+Afc - Xfc) j = P''5s.(Xfc+Aft-Xfe)=0 VO<fc<r - 1 • 

If Sq = then s ■ Xj. = VA;, the Kronecker 6 gives 1 and from Eq. ( pTD we obtain 
|^P(so, Sa = 0) = (P'' — l)5so=o- If on the other hand 7^ then (s ■ X.k)k has maximum 
period and the Kronecker 6 vanishes unless Ak = 0. In this case Eq. (|ll]) yields P^ / [P^ — 

1) - K=o- n 

The choice of parameters is determined by avoiding small s with = 0. For practical 
purposes it is more convenient to replace the condition = (s-Xfc)o<fc<r = by the equiv- 
alent requirement = (s • Xi_fc)i<fc<,,. J2]=k ^j^j-k = mod P , k = 1,2, . . ., min(r, n). 
li n < r the only solution is s = mod P. For n > k one has the problem of finding the 
smallest lattice vector of an n-dimensional lattice. The unit cell of this lattice has the 
volume P^ (cf. Sec. |^). Thus for proper parameters the quality of the generator is 

Nx = P'--l, Q^ = Q2 = ... = Qr = V2/(l- P~') , Qn>r ~ P'/""' • (79) 

This is the first generator which has Qn > 1 for 1 < n < r > 1. For 2 < n < r the 
generator has higher iV„ = max(A^xQn; Nx) = P^ — 1 but lower = max(A^xQn; M) = 
P than the multiplicative generator modP'^ (Ex. 4.1 2, with ~ Mn ~ p**/"). 

In particular if one needs the full periods this generator may be recommended. For 
prime numbers of the form P = 2^ ± 1 the generator has good performance, too. If the 
prime factors of P'' — 1 are known it is no problem to find multipliers which lead to a full 
period. A short discussion of the choice of parameters for large P and r is given in the 
next section and an implementation is presented in Ex. 5.2 2. 



5 Choice of parameters 

We start with a discussion of the mixed multiplicative generator (the multiplicative gen- 
erator is analogous). For practical purposes we can restrict ourselves to M = 2^^ and 
a = 5 mod 8. We set c = 1 which is equivalent to any other odd c and assume n>2 since 
the case n = 1 depends only on b which is 4 for a = 5 mod 8. 
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From Eq. ( |60D we obtain, as long as hsa,M/b < M, that g vanishes unless Sa,M/b\so ^ 
and therefore (A) Qn{sQ,s) = ^1 + s^/s^ M/bl^ ^ -^Z^- However if Sa,M/b = M/b we get 
l^fp = M for (B) So = Sa = OmodM. This leads to Qn = |s|/M which for some s is much 
smaller than 1/4. 

So Eq. (B) is more important. We solve it for si yielding si = kM — as2 — o^s^ — 
Sn depending on the free integer constants k, S2, S3, . . . , s„ which give rise to 



— a 



n—l 



an n-dimensional lattice (cf. 0). The lattice is given by an n by n matrix A according 
to s = A ■ {k, S2, . ■ ., Sn)"^ and we read off 



A 



( M -a -a? 
1 

1 

V 



1 ) 



( M -a 

1 -a 

1 
V 



(80) 



where zeros have been omitted and both matrices define the same lattice since they differ 
only by SL{n, Z) lattice transformations. 

We denote the length of the smallest non- vanishing lattice vector by z/„. Since the 
quality of the random numbers is determined by I'n = MQn we search for an a which 
large i>n. Most important are small n, in particular the pair correlation n = 2. In the 
best case the lattice has a cubic unit-cell and Un is determined by the dimension of the 
lattice and the volume of the unit-cell. Since the volume is given by the determinant of 
A we get as an approximate upper bound z/„ £ M^/". The calculation of z/„ is a standard 
problem in mathematics for which efficient algorithms exist [|T0[] . 

To simplify the search for reasonable multipliers it is useful to have also a lower bound 
for Un- Due to the specific form of A it is easy to see that z/„ has to be larger than the 



smallest ratio > 1 between two elements of then set {1, a, a mod M, 



mod M,M}. 



since ^ M + M^/^ + iM^/^ 



If we take e.g. a « 

we find 1^2 ^ M^/^ which is identical with the upper bound. 
Similarly we obtain ^ M^/^ if we take a ~ M^/^ or a ~ M^/^. However this 
is not compatible with a ^ M^/^ and we only get z/2 ^ M^/^. On the other hand we 
can take a ^ M^/^ + ^M^/^ which differs little from M^/'^. Therefore z/2 ^ M^/^ ^nd 

modM we have z/3 ^ M^/^. Generally, with 
1 A#i/2'=-i (^^Yie plus signs may as well be replaced by minus 
signs) we get ^ M^^^" ^ as long as k > n and M^/^" ^ ^ 1. Note that M^/^" ^ is only 
a lower bound for z/„. In the generic case z/„ will be close to M^/" (cf. Ex. 5.1 1). 

Obviously increases with M. For all practical purposes the magnitude of M is 
only limited by the performance of the generator. In practice one has to split M into 
groups of digits (16 or 32 bit) that can be treated on a computer. The multiplication by 
a performs best if the pre-factors 1/j are omitted. This should be done even though for 
a _^ + . . . + the lower bounds for z/„ decrease, i/„ Z M^/"^"'^ /{n - 1)!. 

Note that the number of digits of M is much more important for randomness than the 
fine-tuning of a. 

Finally, we have to add not too small a constant oq = 5 mod 8 (16 or 32 bit) to the sum 
of powers of M. This constant can be fixed by explicit calculation of the or by looking 
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at (A) from the beginning of this section which imphes that should have large |s| for 
all 16 < m = Sa^lM^^^"'' ■ A suitable choice is e.g. = 3 580 621 541 = 62 181 mod2i^ 
With this value of oq we find |s| m^^^ for n = 2, 3. 

We summarize the result for the parameters of the mixed multiplicative generator: 



M -- 
ao = 
Nx 



C = 1 

5 mod 8 , ao 



, a = 
72/4 



e.g. 

Qn 



ao 



*-'^do _|_ _|_ 22(io 

3 580 621 541mod2'^« , leads to 



22'=-ido ^ 2''^ ''JO ^ _ _ _ ^ 2"'*« + ao , with 



1 



(81) 



Now we turn to the improved generator of Sec. [1.4|. The discussion of the generator 



of Sec. is analogous. 

The Fourier transform of the generator is given by Eq. (|73|). We set n > 2 since 
independently of the parameters Qi = 1. Further on, we fix an m|M and find that = 
m if and only if (C) Sa = km, k odd if m < M, and (D) sq + ^ YJj=2i.^^~^ ~ ^)^3 = 
(We neglect here that |gp may even be 2M for m = M, cf Eq. (^5]).) Eq. (C) can be solved 
for si and Eq. (D) for so depending on the integer parameters fc, £, S2, . . . , s„. This gives 
rise to an {n + l)-dimensional lattice (for m < M we actually obtain an affine sub- lattice 
since k has to be odd) determined by the matrix B via [sq, s) = B ■ {£, k, S2, ■ ■ ., s„)^, 



B 



( m — c 
m —a 
1 



-c(a"-2 + ... + i) / 



m — c — c 
m —a 

1 —a 



\ 



m 



m 



-c 
-a 



a 

-a — 1 
1 



-a 



a + a 
2 - a-1 



/ V 



■ . —a 
1 



(82) 



+ a 



\ 



1 



^3) 



where again zeros have been omitted. 

B describes an [n + l)-dimensional lattice which has a unit-cell with volume m^. However 
this does not imply that the smallest lattice vector Un has length of about m^/*^"^^^ We 
see from (|83D that there exists an {n — l)-dimensional sub- lattice with sq = and S2 = 
—si — S3 — ... — s„ (delete the first and the third row and column in (^3])). The unit-cell of 
the sub-lattice has volume m and z/„ ^ 77i^/("-~i) which, for n > 4, is smaller than m^/^""*"^-'. 
The smallest lattice vector for n > 4 will have the form (0, Si, — Si — S3 — . . . — s„, S3, . . ., s„) 



with the length {sl + sl + . 



-|- s^ -|- (si -|- S3 -|- . . . -|- Sn)"^)^^"^- Since this is of about the same 
magnitude as {s'f + s'^ + . . . + s^)^/^ we may simply omit S2 and reduce the problem to 
the [n — 1) dimensions given by (si, S3, . . ., s„). Geometrically this means that the lattice 
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corresponding to B for n > A never has an approximately cubic unit-cell. Note moreover 
that the sub-lattice is independent of c which means that c can not be fixed by looking 
at the n-tupel distributions for n > 4. 

The smallest value of Qn = Vn/^ is obtained for m = M which is thus the most 
important case. For m = M we are not restricted to odd k. The situation is similar to the 
(n — l)-dimensional case of the mixed multiplicative generator, Eq. (|80|), with —a^ replaced 
by a^+a^~^ + . . .+a. This allows us to use a = M'^/'^ +M^I^ + . . .+M^/^'''^ +ao again. Since 
1 < a ^ Ml/2 ^ a^modM ^ 2M^/^ < . . . < a'^-^modM ^ (n - 2)!Mi-2''" < M 
we have + a^~^ + . . . + a ^ a-^ modM. The minus sign is irrelevant, thus we can 
copy the corresponding lower bounds from the mixed multiplicative generator: ^ 
M"'^/^" /(n — 2)! for + 1 > n > 4. The constant oq is given by the case m < M as will 
be discussed below. 

The constant c can be fixed by the case n = 2. We have to meet two equations (E) 
Si + as2 = OmodM and (F) Sq + CS2 = OmodM to get Ig]"^ = M. Both equations are 
solved by e.g. Sq = — c, si = —a ^ — M^/^, S2 = 1 with |(so,s)| ^ (c^-FM)^/^. In order to 
reach the theoretical limit 1^2 ~ M^/^ one needs c ^ M^/^. So, the simplest ansatz for c is 
c = + 1 for A > 2/3, e N. On the other hand, if Si = M^-^s[, S2 = M^'^s'^ then 
s'-^ + as'2 = mod M'^ has a solution with |s'| £ M'^/^. Since (F) is solved by Sq = —M^'^s'^ 
we find |(so, s)| ^ \sq\ ^ M^'^M'^^'^ = M^'^/"^ . To allow for the maximum value M^/^ one 
needs A < 2/3. In general, c should not have more successive zero digits than M^/^ has. 
The simplest reasonable choice is therefore c = M'^/^ + 1. We can generalize this slightly 
to c = (2'^i + l)co, where cq is a 16 or 32 bit number and 2'^^ < M'^/^ < cq2'^^ . This choice 
of c leads to best performance among all reasonable c. We will see in Ex. 5.1 1 that it 
actually gives 1^2 ~ M'^/^ and 1/3 ~ M^/^. As a lower bound for z/2, z/3 one has only the 
values M-'^/^, M^/'^ that are obtained from Eqs. (E), (C) alone. 

Now we determine Cq and by looking at Sa^M = m < M. The case m < M is 
more important than for the mixed multiplicative generator since Qn is not limited by 
1/4. To some extent the smaller Qn for m < M is compensated by the fact that for 
small m there are more points with Sa = odd ■ m. We use Cq = 3 580 621 541 as for the 
mixed multiplicative generator and find with cq = 3 370 134 727 = 11 463 mod 2^^ that 
Q2{so, s) ^ m2/3-i and Qsiso, s) ^ m^/^-^ if Sa,M = m. 

We summarize the result for the generator of Sec. |4.4|: 



M = 2^''^", a = 22'"^'^" + 2^'-""'^ + ... + 22^0 + ao, c=( 2-t(2'=+V3)do + A ^^^^ 



ao = 5 mod 8 , ao ^ 2"^" , cq odd , co > 2'^^^"^° , 

e.g. ao = 3 580 621 541mod2'^« , co = 3 370 134 727 mod 2'^° leads to 

Nx = 2^'>\ Qi = l, Q2^M^/^-\ g„>3 ^ mV(«-i)-i . (84) 

Finally we give a short discussion of the multiply recursive generator of Sec. ^]5| (cf. 
Ex. 5.1 2). 

P should not be taken too small to provide enough digits for the random numbers. To 
optimize the performance one should use a prime number of the form P = 2'^ ± 1, e.g. 
P = 2'^^ — 1. Moreover we set a^-i = 1, ar-2 = . . . = ai = 0. 
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The most severe problem is to find the prime factors of P*" — 1. To this end it is useful 
to take r = 2^ since in this case P^^ — 1 = (P^* ^ + 1) ■ . . . ■ (P + 1) ■ (P — 1) and one is 

ofc — 1 

basically left with the problem to determine the prime factors of P +1. 

Afterwards it is easy to find an ao € that makes the polynomial P(A) = A'' — A*""^ — 
ao primitive over Zpr. Since 



/ 



\ 



\ Xk-r+1 J 







with X = 



ttr-l Or-2 
1 



V 



ai ao 



1 



(85) 



a necessary and sufficient condition for a maximum period is X^^'"^) = 11 modP and 
X(^''~^yp ^ ImodP for all prime factors p of P^ — 1. High powers of X are easily 
computed. If N = E&i2\ h e {0, 1} then X^ = U{i;b.=i} and X^' = {X'^'''f. 

Now one has to check the n-tupel distributions for n > r. We found (Eq. (|7|)) 
that = P'^ — 1 if and only if sq = and = 0. The latter equation is equivalent 
to = YJj=k^jXj-k = ^fc-P, k = l,2,...,r, ik £ ^ (cf. Sec. |0|) and gives rise to an 



n- dimensional lattice determined by C via s = C ■ {i 



1) 



i) , C — Ci ■ ■ ■ Cr 



( 1 



\ 



P -Xi 

1 



-X 



n—k 



V 



( p 



c ~ Co = 



J 



-ao 



P -1 



V 



1, ar-2 



-ao 



1 ; 

(86) 

= Oi = 0. The 
i>2r = y which 

. Since the 



(after some lattice transformations) if r < n < 2r and a^-i — -l, 
determinant of Co is P*", however the symmetry of Co leads to z/^+i = . . 

( P 

is given by the shortest lattice vector of the 2 by 3 matrix 

\-ao -\ I ) 

second and third row are identical (up to a minus sign) the problem is analogous to the 
calculation of z/2 in the case of the mixed multiplicative generator. We obtain v £ 2^/^P^/^ 
with ao ^ 2i/^pi/2 ^ 55^99 for P = 2^1 - 1. 

We summarize the result for the generator of Sec. 



1, prime. 



1, ar-2 



oi = 0, ao ^ 2^/^P^/^ with 



X^P^-^) = ImodP and X^^''"^)/^ ^ ImodP Vp| (P*^ - I) , p prime, leads to 



X- 



X 



(P" 



Qi 



r+1 



Q2r ^ 2"^P^'^-' . 



^7) 



Notice that the effort for calculating random numbers does not increase with r. 

Let us finally mention that the quality of the n-tupel of the (non-successive) random 



numbers X/, X-, 



k„+e 



deteriorates to Qn ~ P^^ '^^^"^ if there exist d > r—n values 
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of j G {—1, • . ., —r} with Xj = (see the remark at the end of Sec. p.2| ). In particular 
if Xk-i = . . . = Xk-r+i = the pair (Xo,Xk) has quahty of less than pi/^-'' because 
aXi = hXk+i V£ if a = hXk mod P. From Eq. (^) we see immediately that this happens 
for multiples oik = {P"^ — 1)/ {P — 1) (notice the equidistant zeros in Fig. 121). This makes 
it not desirable to use more than (P** — 1)/(P — 1) multiply recursive random numbers. 
Example 5.1. 

1. We set M = 2^^^ = 2^'-^^, a = 2^28 + + 2^2 + 62 181 and in case of the generator 
of Sec. fl.4| c = (2^^° + 1) ■ 11 463. In the following table we compare the mixed 
multiplicative generator with the generator of Sec. |4.4| . The results can easily be 
obtained with a computer algebra program and Eq. ([75|). 





Xk+i — aXk - 


h 1 Eq. (|1|) 


Xk+i = aXk + cmk{k/2) 


Eq. (H) 




0.99414 


0.99414 


1.00000 


1.00000 


a2 


0.50000 


0.50000 


0.65658 


0.66667 


as 


0.33203 


0.33333 


0.49783 


0.50000 


0:4 


0.24859 


0.25000 


0.33436 


0.33333 


0:5 


0.19721 


0.20000 


0.24636 


0.25000 


ae 


0.16335 


0.16667 


0.19882 


0.20000 



(88) 

We see a good agreement of the quality parameters with the approximate upper 
bounds. This means that our choice of parameters is satisfactory. Moreover the 
table confirms that the quality parameter of the generator of Sec. [4.4| lies above the 
quality of the mixed multiplicative generator. 

Finally we present an implementation of the generator in Pascal. We group the 
digits of Xk to 16 blocks of 16 digits X[l] , . . . , X[16] starting from the highest 
digits. 

unit randoml ; 
interface 

const n=16; n0=(n+2) div 3; a0=62181; c0=11463; 
var X : array [1 . .n] of longint; 
procedure nextrandom; 
implementation 

var even: boolean; i:word; c: longint; 
procedure nextrandom; 
var j ,k:word; 
begin 

if even then inc(c,cO); even:=not even; 
for j:=l to n do begin 
X[j] :=X[j]*aO; 

k: =2; while j+k<=n do begin inc(X[j] ,X[j+k] ) ;k:=k shl 1 end end; 
inc(X[n-l] ,X[n] shr 16); X[n]: = (X[n] and $FFFF)+c; 
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inc(X[nO-l] ,X[nO] shr 16); X[nO] :=(X[nO] and $FFFF)+c; 
for j:=n downto 2 do begin 

inc(X[j-l] ,X[j] shr 16); X[j] :=X[j] and $FFFF end; 
X[l] :=X[1] and $FFFF 
end; 

begin for i:=l to n do X[i] :=0; c:=0; even:=true end. 

The corresponding mixed multiphcative generator is obtained by omitting or chang- 
ing the hnes containing c. On a lOOMHz Pentium computer this (not optimized) 
program produces 19 563 random numbers per second whereas 20 938 mixed mul- 
tiphcative random numbers can be produced. A loss of speed of about 6.6% seems 
us worth the gain of better random numbers. Note that the number c suffers an 
overflow every about 750 000th random number. This does not affect randomness 
and it is not worth the effort to correct this flaw. 

2. We set P = 2^^-l,r = 8 which leads to P'^ - 1 = 2^4 . 3^ ■ 5 ■ 7 ■ 11 ■ 17 • 31 ■ 41 ■ 151 • 
331 • 733 ■ 1709 • 21529 ■ 368140581013 • 708651694622727115232673724657. Moreover 
we take a^-i = 1, 0^-2 = . . . = Oi = 0, Oq = 60 045 yielding 

Xo = 1, X_i = . . . = X_7 = , Xk+i = Xfc + 60 045Xfc_7 mod 2^^ - 1 , (89) 

Nx = P«-l « 2'^', Q,^...^Qs^ (p8y.00202-l^ Q^^ ^Q^^^ (p8^ 0.06368-1 _ 

The following program gives on a lOOMHz Pentium 74 473 random numbers (X [k] ) 
per second. 

unit randoin2; 
interface 

var X : array [0 . .7] of longint ; k : integer ; 
procedure nextrandom; 
imp 1 ement at i on 
const a0=60045; 

var i: integer; xO, xl,x2: longint ; 
procedure nextrandom; 

begin 

xO:=X[(k+l) and 7] ; 

x2:=(x0 and $FFFF)*aO; xl:=(xO shr 16)*a0+(x2 shr 16); 
x2:=(x2 and $FFFF)+(xl shr 15)+((xl and $7FFF) shl 16); 
if (x2 shr 31)=1 then x2:=(x2 xor $80000000)+l ; 

inc(x2,X[k]) ; 

while (x2 shr 31)=1 do x2:=(x2 xor $80000000)+!; 
k:=(k+l) and 7; 

if x2=$7FFFFFFF then X[k] :=0 else X[k] :=x2 

end; 

begin k:=0; X[0] :=1; for i:=l to 7 do X[i] :=0 end. 
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6 Results and outlook 

We have generalized the spectral test. As the new feature we analyze the sequence of 
random numbers (I in the figures) not only the distribution of ra-tupels (II in the figures). 

We saw that the mixed multiplicative generator did not pass the test with an ideal 
result. We were able to construct an improved generator which has the recursion formula 

Xo = 0, Xfc+i = aXfc + cint(fc/2) mod2'^. (90) 

For the choice of the parameters a, c, d we made suggestions in Eq. (|84D . This generator 
(or the multiply recursive generator given in Eq. (^)) seems us to be the best choice in 
quality and performance. An implementation of a generator of this type with modulus 
2^ _ 2256 ^ ]^g77 presented in Ex. 5.1 1. The calculation of random numbers is fast 
even though the modulus is that large. We think that for all practical purposes pseudo 
random numbers generated with this generator can not be distinguished from a true 
random sequence. 

We were able to analyze this generator and several others by virtue of a remarkable 
formula on the harmonic analysis of multiplicative rings of remainder class rings (cf. Thm. 
3.4 and Eq. (|47|)). The choice of parameters was discussed in Sec. |. 

For practical purposes there is essentially no need for further improvements. From a 
purely mathematical point of view however there are lots of open questions. 

First of all one could be interested in the cases where Eq. ( |17| ) does not provide the 
result (take e.g. M = P prime and a no primitive element of Zp, cf. Ex. 4.1 5). In these 
cases N\g\'^{so, Si) is given as zero of the polynomial 

PsoiY)^ n {Y-N\g\'iso,s,)) . (91) 

For multiplicative generators Ps^{Y) = — MNY^^~'^ + . . .. Numerical calculations 
show that has integer coefficients. We were not able to prove this for sq 7^ nor to 
analytically determine the coefficients for non-trivial examples. 

Further on, generators involving polynomials may be of interest. An example with a 
quadratic polynomial was presented in Ex. 4.6 4. 

Finally we would be interested in multiply recursive generators. Those generators are 
given by a matrix-valued multiplier. The simplest example with a prime number modulus 
was presented in Sec. [4.5| . In this section we saw that multiply recursive generators are 
also the best candidates for being even more efficient than the generator given in (|90D . 
In this connection multiply recursive generators with power of two modulus may be of 
special interest. 
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Figures 

Some graphs of random number generators are presented to give a visual impression of 
what the generator looks like. There are two possibilities to draw a two-dimensional 
plot: first (I), to plot the k-th random number Xk over k and second (II), to plot X^+i 
over Xk presenting the pair correlation. The third part of the figures give the absolute 
of the Fourier transform of I. \g\'^{sQ,si) is a measure for the correlations along a line 
perpendicular to (sq, si) in I (cf. Eq. ([16|) ). For ideal generators \g\'^ should be < 1 and 
Figs. I and II should look like first rain drops on a dry road. 
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Fig. 3: Xo = 1, Xk+i 
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Fig. 4: Xo = 0, Xk+i 
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Fig. 7: Xo = 0, Xfe+i 
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Fig. 8: Xo = 0, Xk+i = (37 + 1024)^^ + 1 mod 1024^ 
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Fig. 9: Xo = 0, Xk+i 
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Fig. 10: Xo = 0, Xk+i = 30Xfe + 25A;mod 11^ 
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Fig. 12: Xo = 1, X_i = 0, Xk+i ^Xk + 7Xfc_imod31 
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